Filtering a distribution simultaneously in real and Fourier space 
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We present a method to filter a distribution so that it is confined within a sphere of given radius 
r c and, simultaneously, whose Fourier transform is optimally confined within a sphere of radius 
fc c . Our procedure may have several applications in the field of electronic structure methods, like 
the generation of optimized pseudopotentials and localized pseudocore charge distributions. As an 
example, we describe a particular application within the SIESTA method for density functional 
calculations, in removing the spurious rippling of the energy surface generated by the integrations 
in a real space grid. 

PACS numbers: 71.15.-m,31.15.-p,31.15.Pf 



I. INTRODUCCTION 



It is well known that the mean quadratic widths, in 
real and Fourier space, Ar 2 and Afc 2 , of a distribution in 
a space of nr> dimensions, obey the uncertainty relation 
ArAfc > no- The equal sign applies to a spherically sym- 
metric gaussian distribution, which is therefore optimally 
confined in phase space in the least squares sense. In 
many practical cases, however, we are interested in distri- 
butions that are strictly confined within a sphere of given 
radius (i. e. defined to be strictly zero outside that sphere) 
and, simultaneously, optimally confined within another 
sphere in Fourier space. This occurs when, to be com- 
putationally efficient, we use distributions defined only 
within a finite sphere and we must limit also their Fourier 
transforms to a finite number of plane waves. In order to 
calculate those Fourier components, we frequently must 
perform a discrete Fourier transform using a finite num- 
ber of grid points, and we want to avoid as much as possi- 
ble the resulting aliasing effects 1 . Such a situation occurs, 
for example, in the efficient computation of Ewald sums, 
and in the particle- mesh method 2 . Within the field of 
electronic structure calculations, this problem occurs in 
the real-space formulation 3 of Klcinman-Bylander pseu- 
dopotentials 4 , and of pseudocore charge distributions of 
ultrasoft pseudopotentials 5 . 

In the specific case of the SIESTA density functional 
method 6, , this problem arises in the evaluation, using a 
real-space grid, of matrix elements involving strictly lo- 
calized basis orbitals and neutral-atom potentials. Those 
integrals generate an artificial rippling of the total en- 
ergy, as a function of the atomic positions relative to 
the grid points (the so-called eggbox effect), which com- 
plicates considerably the relaxation of the geometry and 
the evaluation of phonon frequencies by finite differences. 
In other grid-based methods 8 , this problem is generally 
solved by filtering the atomic pseudopotentials 9 , typi- 
cally by multiplying them by an ad-hoc filter function 
in Fourier space 1 . Here we present a new method for 
optimal filtering and its application to solve the eggbox 
problem in SIESTA. 



II. OPTIMIZED FILTERING METHOD 

We will study only the specific case of three dimen- 
sions, but the extension to one or two dimensions is ob- 
vious. Consider an initial distribution of the form 



F(r) = l F W y "( f ) ifr ^ • 
^ ' 1 otherwise 



(1) 



where we are using the same symbol F for F(r) and its 
radial part F(r), since it does not lead to any confusion. 
Y" l (r) is a real spherical harmonic. The Fourier trans- 
form of F(r) is 



G(k) 



(2tt) 3 / 2 



dr 3 



r F(r) = G(k)Yr(k) (2) 



where we have introduced the factor v to make G(k) real, 
and 



G(k) = 



4tt 



f'I'c 



(2tt) 3 / 2 



dr r 2 ji(kr)F(r) 



(3) 



where ji (x) is a spherical Bessel function. 

In general G(k) will be nonzero for any value of k. If we 
want to filter it out for k > k c , the most straightforward 
procedure is to multiply it by a step function and then 
to perform the inverse Fourier transform: 



F(r) 



-in 



(2tt) 3 / 2 



dk k 2 ji{kr)G{k). 



(4) 



The new F{r) will no longer be strictly zero for r > r c 
but we may suppress those components and iterate the 
procedure. As a result, only the most confined compo- 
nents, in real and reciprocal space, will survive. 

To annalize more rigorously the decomposition of F{r) 
into more and less confined components, let us define 
x = r/r c , y = k/k c , f(x) = xF(x r c ), g(y) = yG(yk c ), 
k = k c r Cl and K{x,y) = nxy ji(nxy). Then, 

substituting in (3) and (4) , one iteration of the filtering 
procedure is given by 



/(a 



dx'K 2 (x,x')f(x') 



(5) 



2 



where 



if 2 (x,x 



Jo 



dyK(x,y)K(y,x') 



(6) 



If /(x) were already a perfectly confined function in 
both real and reciprocal space, it would not be affected 
by the filtering procedure (5), i. e. it would be an eigen- 
function of the filtering kernel K 2 with eigenvalue one. 
In practice, the uncertainty principle forbids simultane- 
ous perfect confinement in real and Fourier space, and 
the filtered /(x) will unavoidably 'leak' somewhat out- 
side x > I and its norm within x < 1 will no longer be 
one. In fact, if <fi(x) is an eigenfunction of K 2 , with norm 
equal to one within x < 1, its eigenvalue A 2 gives directly 
its norm after filtering, since the effect of filtering is just 
a multiplication by A 2 : 

4>(x)<- [ dx'K 2 (x,x')<f>(x') = \ 2 <P{x) (7) 
Jo 

Thus, we may perform an efficient filtering, without the 
need of iteration, by expanding the original function in 
terms of the complete basis of cigenfunctions of K 2 , keep- 
ing only those with eigenvalues sufficiently close to one. 
Since, it is clear that the eigenfunctions 4>i{x) of K{x, y), 
with eigenvalues Aj, are also eigenfunctions of K 2 (x,x r ), 
with eigenvalues A 2 , we may work with the simpler eigen- 
value problem 



/' 

Jo 



dyK(x,y)4>i{y) = A^x) 



(8) 



Notice that, since K (x, y) is the Fourier-transform kernel, 
the eigenfunctions (pi (x) have the same shape in real and 
reciprocal space. This is not true in general for the fil- 
tered function /(x), which is a combination of eigenfunc- 
tions with eigenvalues A^ close to either +1 or -I, which 
either change sign or not when Fourier transformed. 

In order to solve (8), it is convenient to expand K (x, y) 
and <pi(x) in a basis of functions in the interval [0,1]. 
The simplest basis is that of powers of x. From the 
Taylor expansion of ji{x) at x = we find K(x,y) ~ 

En=o K nX 2n+l+1 y 2n+l+ \ where 



K n 



2k {-\) n K 2n+l+l 



(9) 



7T (2n)!!(2n + 2Z + l)!! 
Then making 0,(x) = J2n=o 4>mX 2n+l+1 , Eq. (8) becomes 



N 

E 

m=0 



2n + 2m + 21 + 3 



A»<; 



(10) 



In practice, we have found numerically more accurate, 
stable, and efficient (requiring a lower N) to expand 
K(x,x') in orthonormal Legendre polynomials 1 P n (x) 
in the interval < x < 1. Taking into account the parity 
l p = mod (1,2) of ji(x): 

N 

K{x,y)~ ^2 K nm P 2n -i p -i{x)P 2m -i p -i{y). (11) 

n,m— 1 



The kernel coefficients K nm may be calculated by inte- 
gration in a Gauss-Legendre 1 set of points x Q and weights 



K nm = J J dx dyK(x,y)P 2n -i p -i{x)P 2m - 



y(y) 



N- 



= ^ w aWpK(x a ,yp)P2n-l p -l{x a )P2 m -l p -l{ty$) 

The required number TV of polynomials is deter- 
mined by the convergence of the expansion xji(x) ~ 
J2 n =i 3lnP2n-l p -i{x) in the interval < x < k. Fig- 
ure 1 shows the number of polynomials N required to 
obtain a given error in the expansion, as a function of k, 
for I = 0. The Z-dependence of the error is very small 
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FIG. 1: Number of polynomials required to obtain a root 
mean square error of the expansion of xjo(x) in the interval 
< x < k c r c . jo(x) is a spherical Bessel function with I = 0. 

and, as a rule of thumb, we use N — int(10 + 0.65/t). 

Figure 2 plots the first eigenfunctions </>j(x) of the filter 
kernel K 2 (x, y) for a typical value of k, and figure 3 shows 
all the eigenvalues A 2 up to N. It may be seen that there 
is a rapid transition between the eigenvalues which are 
very close to 1 and those close to 0. It is then straight- 
forward to select the M eigenfunctions whose eigenvalues 
are above some threshold, say A 2 > 0.99, for the expan- 
sion of the filtered function: 



M 



fix) <- ^2fi<t>i{x) 



i=i 



N- 



fi = J! w a<f>i{Xa)f(x a )- 



(13) 



(14) 



a=l 



Fig. 4 shows, as an example, the unfiltered and filtered 
oxygen 2p pseudo atomic orbital, generated as proposed 
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FIG. 2: First few eigenfunctions (with highest eigenvalues) 
of the filter kernel K 2 (x,x') for k = k c r c = 25. Divided by 
r, they give the radial part of the distributions, with angular 
momentum Z, that are most localized in a real-space sphere 
of radius r c and simultaneously in a reciprocal-space sphere 
of radius k c . 
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FIG. 3: Eigenvalues of the filter kernel K 2 (x,x') for k — 25 
and I = (circles and full line), I — 1 (squares and dashed 
line), and I — 2 (diamonds and dotted line). 
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FIG. 4: Filtered (dashed lines) and unfiltered (full lines) oxy- 
gen 2p pseudo atomic orbital, generated with a strict cutoff 
r c — 3.94 bohr as proposed by Sankey and Niklewski 10 . It 
was filtered with a cutoff k = 25, which corresponds to a 
plane wave cutoff k c = 6.35 bohr -1 , or k\ = 40 Ryd. Upper 
panel: real space shape. Lower panel: tails of their Fourier 
transform. 



by Sankey and Niklewski 6 ' 10 with a Troullier-Martins 
pseudopotential 11 . To enhance the filtering effect, we 
have used a very small filter cutoff. Still, it may be seen 
that the Fourier components above the cutoff are very ef- 
ficiently suppressed, although this is achieved (with this 
small cutoff) at the expense of a substantial change in its 
shape. 

Finally, the most confined eigenfunctions, <f>i(x), for 
each angular momentum I, may be used to generate a 
localized distribution with given multipole moments, as 
required in the ultrasoft pseudopotential 5 and projector 
augmented waves 12 methods, among others problems in 
computational physics 2 . They may be used also as a 
basis of localized orbitals, for the expansion of the elec- 
tron wavefunctions 13 , which is asymptotically complete, 
within the confining spheres, as the filter cutoff increases. 
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III. APPLICATION WITHIN SIESTA 



There are three contributions to the eggbox effect in 
SIESTA (an artificial rippling of the total energy sur- 
face as a function of the positions of the atoms relative 
to the integration grid points): i) the so-called neutral- 
atom potential 6 Vna( v ) given by the local part of the 
atomic pseudopotentials minus the Hartree potential of 
the free-atom electron densities; ii) the exchange and cor- 
relation potential V xc (r), given by the electron valence 
density /o(r), which in turn is given by a sum of products 
of atomic basis orbitals f^v). These two contributions 
are frequently comparable in magnitude; in) the nonlocal 
core correction (NLCC) to V xc (y), given by a pseudocorc 
electron density Pnlcc{ v ) added to p(r). This added 
density is generally very large and localized and, when 
the NLCC is present, it normally dominates the eggbox 
effect. Finally, the Hartree energy, given by the self- 
interaction of p(r), also contributes to the eggbox but, 
since the Hartree potential is much smoother than the 
density this contribution is always negligible compared 
to the other ones. 



Thus, in order to cut drastically the eggbox effect, we 
must filter PNLCc( r ), Vna(t), and ^(r). The first two 
may be filtered with the plane wave cutoff k c of the real- 
space integration grid used to calculate the matrix ele- 
ments of Vna(?) and V xc (r). The filtering cutoff required 
for <p^(r) is somewhat less clear, because we need to treat 
products of two <^'s in the integration grid, not just the 
y?'s themselves. In principle, the plane wave cutoff of a 
product of two functions is twice that of the functions 
themselves, what would suggest that ipy,{v) should be fil- 
tered with k c /2. However, a widespread experience with 
plane wave codes has shown that this criterion is too 
strict, and that in practice the effective cutoff for the 
density is typically less than two times that of the wave- 
functions. Therefore, we have checked that making the 
filter cutoff for (^(r) equal to ~ 0.6fc c leads generally to 
the best convergence, as a function of k c . 



Figure 5 shows the eggbox effect of isolated atoms dis- 
placed across the integration mesh. It may be seen that 
the effect is indeed eliminated almost completely by fil- 
tering. Of course, we shall not eliminate the eggbox effect 
at the expense of filtering the pseudopotentials and ba- 
sis functions so much as to change the physical results. 
Figure 6 shows the vibrational frequencies of the water 
molecule, calculated by diagonalizing the dynamical ma- 
trix obtained by finite differences 14 . As the plane wave 
cutoff k c of the integration grid is reduced, Pnlcc{ v ), 
Vna{y) are filtered with that cutoff, and y M (r) is filtered 
with 0.7fc c . It may be seen that much lower cutoffs are re- 
quired, to converge accurate frequencies, with than with- 
out filtering. 
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FIG. 5: a) Total energy of an isolated carbon atom, as it is 
displaced in a large unit cell, using an integration grid with 
a plane wave cutoff = 50 Ryd, whose points are separated 
by Ax — n/k c = 0.44 bohr. b) Magnitude of the eggbox 
effect (peak to peak of total energy) for several isolated atoms 
with hard pseudopotentials or nonlocal core corrections (in 
Pb and Fe), as a function of the plane wave cutoff of the 
integration mesh. Full lines: without filtering. Dashed lines: 
with filtering. 
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FIG. 6: Vibrational frequencies of the water molecule, calcu- 
lated from the hessian matrix, which was obtained by finite 
differences from the forces on the atoms displaced from their 
equilibrium positions. The x axis is the plane wave cutoff 
of the integration grid used in SIESTA. Full lines: without 
filtering. Dashed lines: with filtering. 



IV. CONCLUSIONS 

We have presented a general method to generate dis- 
tributions, with a given angular momentum, which are 
optimally confined within a strict cutoff in both real and 
Fourier space. They can be used by themselves, as to 
produce localized distributions with given multipole mo- 
ments, or as a basis for expanding and filtering an arbi- 
trary initial distribution. As an example, we have shown 
how they can be used to filter the pseudopotentials and 
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basis functions in the density functional method SIESTA, 
thus eliminating the cggbox effect on the total energy, 
due to the calculation of matrix elements in a real space 
integration grid. 



Thus, our basis functions 4>(x) are the normalized dis- 
tributions which are strictly confined to < x < 1 (i. 
e. r < r c ) and whose Fourier transform has the smallest 
norm in y > 1 (k > k c ). 
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Interestingly, an alternative variational principle may 
be demonstrated for a basis of gaussian functions. Thus, 
we maximize the confinement of a normalized distribu- 
tion (j)(r) and its Fourier transform 0(k), in the sense of 
least squared dispersion: 



APPENDIX A: VARIATIONAL PRINCIPLES 

Here we show that the filtering basis functions 4>i(r) 
obey a simple variational principle, and we also present 
an alternative principle for gaussian basis functions. It 
may be easily shown, by a straightforward functional 
derivative, that the eigenvalue equation (7) is equivalent 
to the variational principle 

/ [ dx dx'cj)(x)K 2 (x,x')(j)(x') = max (Al) 
Jo Jo 

subject to the condition of normalization of <f>(x) within 
< x < 1. Now, using that the Fourier transform of 
4>{x) is 



g{y) 



dx 4>{x) K(x, y) 



(A2) 



as well as the definition (6) and the fact that the total 
norm of a function is the same in real and Fourier space: 



f dy [ dx<f>(x)K(x,y) [ dx'K(y,x')<j)(x') 
o Jo Jo 

J d y g 2 (y) = !- d y g 2 (y) = max 

dy 9 2 {y) = min (A3) 



r 



! (k) 



(A4) 



where r c and k c are here scale factors that determine 
the relative confinement in real and Fourier space, rather 
than strict cutoffs. Multiplying by k 2 /2: 



1,2 2 

dv- cY ^' 



2rl 



2 (r) + / dk— (j) 2 (k) = min 



(A5) 



Now, the first term is the potential energy of a quan- 
tum harmonic oscillator with spring constant (k c /r c ) 2 
and wave function <f>(r), and the second term is its ki- 
netic energy. Its well known solutions are gaussians times 
Hermitc polynomials 16 . A similar (but not orthonormal) 
basis, made of gaussians times powers of r, was used by 
Hartwigsen et al 17 to generate compact separable pseu- 
dopotcntials. 
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